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Sunspots are prominent manifestations of solar magnetoconvection and imaging their subsurface 
structure is an outstanding problem of wide physical importance. Travel times of seismic waves that 
propagate through these structures are typically used as inputs to inversions. Despite the presence of 
strongly anisotropic magnetic waveguides, these measurements have always been interpreted in terms 
of changes to isotropic wavespeeds and flow-advection related Doppler shifts. Here, we employ PDE- 
constrained optimization to determine the appropriate parameterization of the structural properties 
of the magnetic interior. Seven different wavespeeds fully characterize helioseismic wave propagation: 
the isotropic sound speed, a Doppler-shifting flow-advection velocity and an anisotropic magnetic 
velocity. The structure of magnetic media is sensed by magnetoacoustic slow and fast modes and 
Alfven waves, each of which propagates at a different wavespeed. We show that even in the case of 
weak magnetic fields, significant errors may be incurred if these anisotropics are not accounted for 
in inversions. Translation invariance is demonstrably lost. These developments render plausible the 
accurate seismic imaging of magnetoconvection in the Sun. 



Sunspots are substantial deviations from the quiet 
Sun, with umbral temperatures dropping by as much 
as 20% from ambient conditions. Numerous questions 
swirl around sunspot physics, such as understanding their 
long-time stability (compared to convective turnover 
timescales) and appreciating their creation, emergence 
and eventual death. The use of helioseismic waves to 
probe the structure of sunspots has a long and contro- 
versial history (for a review, see, e.g., Inversions 
for sunspot sub-surface structure and dynamics (e.g., (H) 
attempt to explain away the observed effects on seis- 
mic waves by an entirely isotropic wavespeed, an ap- 
proximation that has faced subsequent marginalization 
(e.g., owing to the widespread recognition of strong 
anisotropics prevalent in sunspots. Forward modeling of 
wave propagation in sunspots has generated a deeper ap- 
preciation for measurements and fully realistic non-linear 
sunspot evolution calculations ([4]) have proven success- 
ful. However, posing an inverse problem that accounts 
for these anisotropics remains an outstanding problem of 
great relevance, towards whose eventual resolution this 
article takes a significant step. 

Waves are excited stochastically in the Sun due to the 
action of vigorous near-surface convection. The cross 
correlation of wavefield velocities measured at the photo- 
sphere of the Sun (by measuring Doppler shifts of absorp- 
tion lines formed at the photosphere; e.g., the Helioseis- 
mic and Magnetic Imager onboard the Solar Dynamics 
Observatory [51) is eropirically known to be an ergodic 
random process (e.g., [6|). The associated travel time of 
a wave between two points on the solar photosphere is 
estimated by fitting the cross correlation of the plasma 
velocities measured at those points. We introduce a mis- 
fit functional, defined as the L2 norm of the difference 
between observed (rf ) and predicted (rf ) travel times 



along a collection of paths i: ^^{r^ — )^. We then 
pose the following PDE-constrained optimization prob- 
lem 

X = E(< - ^')' - j d^j <hjX-{Ci- S), (1) 

where x is the cost function, S the wave source, uo tem- 
poral frequency, x the spatial coordinate and A(x, cj) a 
vector Lagrange multiplier, the dual to the wave displace- 
ment ^(x, cj). The predicted travel times are linear func- 
tionals of wave displacement ^. The helioseismic wave 
operator C comprises temporally stationary model prop- 
erties, which we attempt to determine. We reproduce it 
here in the temporal- frequency domain (e.g., [3]), 

£| = -cjVC - i^P^i - 2iuopM • V| (2) 
- V(cVV • - V(C • Vp) + gV • (pO 

-^(VxB)x[Vx(CxB)]-^{Vx[Vx(CxB)]}xB, 

where the properties of interest are the density p, sound 
speed c, vector magnetic field B and flows v. Wave 
damping is denoted by T and gravity by g, where these 
are held fixed and not considered to be parameters here, 
the justifications for which may be found in Q. Accelera- 
tion, wave damping, and Doppler shifting by flow advec- 
tion are the first two terms of operator ([2]), the isotropic 
wavespeed term and buoyancy terms form the second 
line and the final two terms are due to the anisotropic 
Lorentz force. Background pressure p is constrained by 
the magneto-hydrostatic (MHS) equilibrium equation 

Vp=-/,g + V- (^BB-^I^, (3) 

where I is the identity tensor. Two scalar thermal param- 
eters, density and sound speed, and two vector quantities. 
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flows and magnetic flelds, are the independent variables 
that we invert for and pressure is constrained by equa- 
tion dSl). 

Given wave operator ([2j), we study in detail the vari- 
ation of the cost function ([1]) to obtain the gradients of 
the travel-time misflt functional with respect to model 
parameters (e.g., 0, B])- These gradients, also known as 
sensitivity kernels or Frechet derivatives, are indicators 
of information in the seismic wavefield and their sensi- 
tivity to relevant model parameters. The computational 
realization of this method (0, Isl) for two-point correla- 
tion function measurements, as in the Sun, requires cal- 
culating a predictive so-called forward wavefleld and an 
adjoint wavefleld that assimilates the misflt. Sensitivity 
kernels for various physical quantities (which form the 
corrector) emerge from a temporal convolution of these 
two wave flelds, allowing us to pose an inverse problem 
of the form 

6x = - I ^ixi^c^c + Kv-(5v + KB-(^B + i^;^lnp, (4) 

where Kc^ Kv, Kb are sensitivity kernels for sound speed, 
flows and magnetic flelds, while K'^ the kernel for den- 
sity, termed as an impedance kernel in geophysics jargon, 
is sensitive to reflectors (e.g., [ol). Equation dH states 
that the measured travel-time shift comprises a sum of 
volume integrals of these perturbations weighted by the 
corresponding flnite- frequency wave sensitivities. Travel 
times are very weakly sensitive to density variations but 
record sharp contrasts in impedance, such as (possibly) 
the horizontal boundary of a sunspot. We do not expect 
to be able to image such reflections since the wavelength 
of waves that we consider are large in comparison to pos- 
sible rapid variations, but retain and compute these ker- 
nels in any case (see Figures 2 and 4 of supplemental 
material) . 

A critical aspect to setting up an inverse problem is 
in appreciating the physical variables to which waves are 
sensitive. It is seen that the variation of the operator (j2j) 
has terms (among others) that contain p8c^ , p^v and (5B 
(see also detailed expressions for kernels in [7]). This im- 
mediately tells us that kernels for sound speed and flows 
are weighted by the density of the model, in contrast to 
kernels for the primitive magnetic fleld. Further, vari- 
ables c and V are forms of wavespeed, which suggests 
the use of Alfven velocity, a = B/v^47rp instead of the 
primitive B fleld. One may conceive of it as a descriptor 
of the anisotropic wave velocity to which waves are di- 
rectly sensitive. Straightforward manipulation allows us 
to rewrite the kernels as follows 

(5B = Airp) = Airp Sa + ^^a/ ^'^'^ ^ (^) 
which together with equation (|4]) gives 

v/4^KB = Ka i^;^i^; + iKa-a, (6) 



thus providing a new expression for variations in the mis- 
flt 

Sx = - [ dyiKcSc^K^-Sv^K^-Ssi^K'pSlnp. (7) 

Jq 

We note that the flrst three terms represent three types 
of wavespeeds, an isotropic sound speed, an advection re- 
lated flow velocity and lastly, an intrinsically anistropic 
velocity. Although not shown here, weighting the mag- 
netic fleld kernels by the square-root of density redis- 
tributes incoherent sensitivity from the upper-most at- 
mospheric layers to the photosphere and shallow interior. 
The transformation for the density kernel in equation (|5j) 
now contains a contribution from the Alfven velocity, and 
could in principle be used to image reflections off sharp 
velocity contrasts. 

The single-scattering flrst-Born approximation cannot 
capture the full scope of wave propagation in strong per- 
turbations such as sunspots (i.e., with respect to the quiet 
Sun; e.g., [iq|)- This implies that inversions for the sub- 
surface structure of sunspots are likely to require an iter- 
ative algorithm, since we have to sequentially reflne the 
predicted travel times, which are nonlinear ly related to 
changes in the model. Thus, in the analysis here, we 
construct a 'sunspot' in MHS equilibrium (Eq. [3]) and 
determine sensitivity kernels relative to this model. 

We introduce a 2-D stream function iIj{x,z) such that 
the magnetic fleld is given by B = {—dzil^^dxip)- Since 
g = (0, —g), Lorentz forces in the x direction are solely 
balanced by the pressure gradient in equation ([3]), i.e., 
dxP = dx{Bl/2 - Bl/2)^dz{BxBz). From this equation 
we calculate the pressure distribution required to support 
this fleld conflguration and then use the z component of 
equation ([3]) to obtain the associated density. Generat- 
ing an MHS state is non trivial since density and pres- 
sure decrease exponentially as a function of height above 
the photosphere; consequently, a large range of choices 
for the fleld conflguration results in negative pressures or 
densities or both. Field conflgurations with strong hori- 
zontal and vertical flelds also require the action of flows 
to maintain force balance, an aspect we do not consider 
here because the complexity of such a model renders difli- 
cult the interpretation of the attendant kernels. We show 
one example fleld conflguration in Figure [H 

A major difliculty in simulating wave propagation 
through strong magnetic flelds is that Alfven speed ||a|| 
becomes extremely large in the atmospheric layers of the 
Sun (due to the exponentially rapidly decreasing den- 
sity), resulting in a very stiff differential equation. Fur- 
ther, wave travel times are very weakly sensitive to the 
dynamics of these layers because the modes are trapped 
below the photosphere. A multiplicative prefactor is in- 
troduced to control the amplitude of the Lorentz force 
terms in ([2j), (e.g., 0, [Tlj). However, this method results 
in a model that is not seismically reciprocal (e.g., 0), a 
central requirement in the formal interpretation of helio- 
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FIG. 1. Magnetic field configuration in our calculations. 
Top panels show Aflven speeds ax = Bx/ V^^rp and az = 
Bz/\^4:7Tp, which are signed quantities. The bottom left panel 
is the ratio of the absolute Alfven speed to the local sound 
speed and is seen to be on the order of 1 at the photosphere. 
The field is relatively weak with the highest Alfven speed 
around 35 km/s and a Wilson depression of 250 km. (for an 
expanded view, see Figure 1 of the supplemental material) 

seismic measurements and the determination of sensitiv- 
ity kernels. Here, in order to maintain seismic reciprocity 
while still saturating the Alfven speed at 40 km/s, we di- 
rectly multiply the magnetic field by a prefactor. While 
this results in a background field configuration that has a 
non-zero divergence, we note that small-amplitude oscil- 
lations about this field are still divergence free. Further, 
in the scheme of linear inversions for magnetic structure, 
the divergence-free nature of the background field is not 
a strict requirement but could be considered a regular- 
ization term. 

We perform linear magneto-hydrodynamic (MHD) 
wave propagation simulations in Cartesian geometry, us- 
ing the pseudo-spectral code SPARC ([li-Iil). Hor- 
izontal derivatives are computed using Fast Fourier 
Transforms, vertical derivatives are estimated on a non- 
uniform grid using compact finite differences ([ill) and 
time-stepping is effected through the repeated applica- 
tion of an optimized Runge-Kutta scheme ([l6|). Ver- 
tical boundaries are lined with absorbent convolutional 
perfectly matched layers ([l3) ^hat are designed to ab- 
sorb MHD waves as well. We implement a phenomeno- 
logical wave damping term along the lines of the recipe 
suggested by [l8[ . Because we restrict ourselves to a 2-D 
field configuration in this problem, Aflven waves are dis- 
allowed and only magneto- acoustic fast and slow waves 
propagate. 

We focus here on the diagnostic ability of the surface 
/ and acoustic pi modes, so chosen because of their sig- 
nificant sensitivity to surface layers. The measurement 
consists of ridge filters applied to isolate these modes. 
The sunspot is assumed to be located at disk center, 
implying that the line-of-sight component is co-aligned 
with the (vertical) z axis. Thus the vertical wavefield 
displacement is used to define the cross correlation mea- 
surement. We show the power spectra and cross cor- 
relations in Figure O We employ the linear travel-time 
definition ( [gI, [19| ) , also used previously by [7] to estimate 
travel time shifts from cross correlations. 
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FIG. 2. Expectation value of the power spectrum of the pi 
and / ridge- filtered measurements (top panels). The limit 
cross correlation C{t) between a point 15 Mm from the left 
of the sunspot center to a point 10 Mm on the right of the 
center is shown for the pi measurement (middle panel). The 
/-mode cross correlation is between the center of the sunspot 
and a point 10 Mm to the right (bottom panel). See Fig- 
ures [3] and [4] also. The positive-time branch is sensitive to 
waves that first arrive at one measurement point and subse- 
quently at the other and vice versa. The loss of translational 
variance implies that the absolute locations of the points mat- 
ter. The dot-dash boxes indicate the measurement windows. 
Travel time shifts of waves are obtained by estimating the 
deviation of the cross correlation from a reference wavelet. 
Mean travel times, defined as the average of the time shifts of 
oppositely traveling waves, are thought to be largely sensitive 
to structure. Difference travel times, defined as the differ- 
ence between the shifts, are considered primarily sensitive to 
symmetry-breaking flows, (for an expanded view, see Figure 
2 of the supplemental material) 



Figure [3] (see also Figures 5 and 6 in the supplemental 
material) displays the sensitivity of the surface /-mode to 
the sunspot. Because we model waves as finite spatial ob- 
jects, their sensitivities extend beyond just the ray path. 
It can be seen that the effect of the spot is significant 
in that the kernels are noticeably asymmetric between 
the point-pair. The time shifts induced by the magnetic 
field are considerable, comparable in magnitude to those 
induced by flow and thermal perturbations. There are 
hints of mode conversion from / to pi in the difference 
kernel for sound speed (top), just below the pixel on the 
right. 

In Figure H] (see also Figures 7 and 8 in the supplemen- 
tal material), we show a set of difference pi-mode kernels 
for a point pair separated by a distance of 25 Mm re- 
spectively. Because the magnetic field is relatively weak 
compared to a sunspot, the acoustic pi mode, whose en- 
ergy is focused in the sub-surface layers, is much less af- 
fected by the field than the / mode. Symmetry is nearly 
completely restored to the pi kernels. 

The Alfven speed kernels for both / and pi modes show 
features of high spatial frequency, and contain signatures 
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of fast and slow magneto- acoustic waves. In the umbral 
regions of the sunspot, waves of high spatial frequency 
are seen to be propagating toward the interior (plausibly 
slow waves). 

Our computations support the view that inversions 
for sunspots, especially when using surface modes, are 
greatly over-simplified if anisotropic wave speeds are 
not taken into account. The realization of this method 
has required a number of theoretical and numerical ad- 
vances, paving the way for seismic imaging of magneto- 
convection in the solar interior. 



FIG. 3. /-mode (surface) wavespeed kernels for a difference 
travel-time measurement between a point-pair 10 Mm apart. 
Kernels sensitive to isotropic sound speed, Alfven speed ax 
and vertical flows Vz are shown. The boundary of the spot, 
marked by the solid black line, is much smaller than the hor- 
izontal wavelength. The horizontal dot-dash line denotes the 
height at which observations are made in the quiet Sun and 
the symbols mark the measurement points. The /-mode is 
seen to be significantly affected by the spot, as seen in the loss 
in symmetry of the kernels. Signatures of magneto-acoustic 
slow and fast modes and hints of conversion to acoustic pi 
may be plausibly discerned upon examination. The integrals 
of the kernels show that the travel times are significantly af- 
fected by the presence of even this relatively weak magnetic 
field, (for an expanded view, see Figure 3 of the supplemental 
material) 
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FIG. 4. pi-mode wavespeed kernels for a difference travel- 
time measurement between a point-pair 25 Mm apart. The 
panels from top to bottom show kernels sensitive to sound 
speed (top), Alfven speeds ax and vertical flows Vz. The 
boundary of the spot, marked by the solid black line, is much 
smaller than the horizontal wavelength. The horizontal dot- 
dash line denotes the height at which observations are made in 
the quiet Sun and the symbols mark the measurement points. 
Plausible signatures of slow modes propagating down into the 
tube may be discerned in the middle panel, (for an expanded 
view, see Figure 4 of the supplemental material) 
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FIG. 1. Figure 1 of the main article, reproduced: Magnetic field configuration in our calculations. Top panels show Aflven 
speeds ax = Bx/^/^iip and = Bz/ V^^rp, which are signed quantities. The bottom left panel is the ratio of the absolute 
Alfven speed to the local sound speed and is seen to be on the order of 1 at the photosphere. The field is relatively weak with 
the highest Alfven speed around 35 km/s and a Wilson depression of 250 km. 



Firstly, we show expanded versions of the figures of the main article. 

Figure [5] displays the sensitivity of the surface /-mode to the sunspot. It can be seen that the effect of the spot 
is significant in that the kernels are noticeably asymmetric between the point-pair. The time shifts induced by the 
magnetic field are considerable, comparable in magnitude those induced by flow and thermal perturbations. There 
are hints of mode conversion from / to pi in the difference kernel for sound speed (top), just below the pixel on the 
right. Plausible signatures of magneto- acoustic slow waves are seen in the difference kernel for and the mean kernel 
for Vz- 

Figure [6] shows /-mode kernels sensitive to density. 

In Figure[7l we show a set of pi-mode kernels for a point pair separated by a distance of 25 Mm respectively. We note 
right away that the acoustic pi is less affected by the spot than the / mode. Advective and Alfven speeds possess similar 
features, i.e., the Vx and ax mean travel-time kernels both show lobes of sensitivity near the observation points whereas 
the Vz and ones show broad sensitivity between the point-pair. In contrast to both the advective and isotropic 
wavespeed kernels, the Alfven speed kernels only show features of high spatial frequency, and contain signatures of 
fast and slow magneto- acoustic waves. In the umbral regions of the sunspot, waves of high spatial frequency are seen 
to be propagating toward the interior (plausibly slow waves) whereas waves at significant inclination to the the z axis 
seem to propagate toward the observation points (plausibly fast waves). 

The sound- and Alfven speed kernels are sensitive at different spatial locations to structure, the former being 
restricted to sub-photospheric layers, the latter showing the greatest sensitivity at atmospheric layers where low 
densities imply very large Alfven speeds. Consequently, a fractional change in Alfven speed at these high layers creates 
travel-time shifts larger than a corresponding change in sound-speed. Figure [8] shows pi-mode kernels sensitive to 
density. 



3 




FIG. 2. Figure 1 of the main article, reproduced: Expectation value of the power spectrum of the pi and / ridge- filtered 
measurements (top panels). The limit cross correlation C{t) between a point 15 Mm from the left of the sunspot center to 
a point 10 Mm on the right of the center is shown for the pi measurement (middle panel). The /-mode cross correlation is 
between the center of the sunspot and a point 10 Mm to the right (bottom panel). See Figures [5] and [3 also. The positive- 
time branch is sensitive to waves that first arrive at one measurement point and subsequently at the other and vice versa. 
The loss of translational variance implies that the absolute locations of the points matter. The dot-dash boxes indicate the 
measurement windows. Travel time shifts of waves are obtained by estimating the deviation of the cross correlation from a 
reference wavelet. Mean travel times, defined as the average of the time shifts of oppositely traveling waves, are thought to 
be largely sensitive to structure. Difference travel times, defined as the difference between the shifts, are considered primarily 
sensitive to symmetry-breaking flows. 
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FIG. 3. Figure 3 of the main article, reproduced: /-mode (surface) wavespeed kernels for a difference travel-time measurement 
between a point-pair separated by a distance of 10 Mm. Kernels sensitive to isotropic sound speed, magnetic (Alfven) speed ax 
and vertical flows Vz are shown. The boundary of the spot, marked by the solid black line, is much smaller than the horizontal 
wavelength. The horizontal dot-dash line denotes the height at which observations are made in the quiet Sun and the symbols 
mark the measurement points. The /-mode is seen to be significantly affected by the spot, as seen in the loss in symmetry 
of the kernels. Signatures of magneto-acoustic slow and fast modes and hints of conversion to acoustic pi may be plausibly 
discerned upon close examination. The integrals of the kernels show that the travel times are significantly affected by the 
presence of even this relatively weak magnetic field. 
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FIG. 4. Figure 4 of the main article, reproduced: pi-mode wavespeed kernels for a difference travel-time measurement between 
a point-pair separated by a distance of 25 Mm. The panels from top to bottom show kernels sensitive to sound speed (top), 
Alfven speeds ax and vertical flows Vz. The boundary of the spot, marked by the solid black line, is much smaller than the 
horizontal wavelength. The horizontal dot-dash line denotes the height at which observations are made in the quiet Sun and the 
symbols mark the measurement points. Plausible signatures of slow modes propagating down into the tube may be discerned 
in the middle panel. 
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FIG. 5. /-mode (surface) wavespeed kernels for a point-pair separated by a distance of 10 Mm. The rows display isotropic, 
advective and anisotropic (Alfven) wavespeeds, in that order. Kernels on the left column are for the mean travel time between 
the point pair whereas the panels on the right column shows kernels for the difference travel time. The boundary of the spot, 
marked by the solid black line, is much smaller than the horizontal wavelength. The horizontal dot-dash line denotes the height 
at which observations are made in the quiet Sun and the symbols mark the measurement points. The /-mode is seen to be 
significantly affected by the spot, as seen in the loss in symmetry of the kernels. Signatures of magneto- acoustic slow and fast 
modes and hints of conversion to acoustic pi may be plausibly discerned upon close examination. 
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FIG. 6. /-mode (surface) density kernels for mean (upper panel) and difference (lower panel) travel-time measurements between 
a point-pair separated by a distance of 10 Mm. The boundary of the spot, marked by the solid black line, is much smaller 
than the horizontal wavelength. The horizontal dot-dash line denotes the height at which observations are made in the quiet 
Sun and the symbols mark the measurement points. The /-mode is seen to be significantly affected by the spot, as seen in the 
loss in symmetry of the kernels. The kernel shows significant complexity, with the difference kernel possessing relatively larger 
sensitivity to density variations within the spot. 
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FIG. 7. pi-mode wavespeed kernels for a point-pair separated by a distance of 25 Mm. The rows display isotropic, advective 
and anisotropic (Alfven) wavespeeds, in that order. The kernels on the left column are for the mean travel time between the 
point pair whereas the panels on the right column shows kernels for the difference travel time. The boundary of the spot, 
marked by the solid black line, is much smaller than the horizontal wavelength. The horizontal dot-dash line denotes the height 
at which observations are made in the quiet Sun and the symbols mark the measurement points. 
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FIG. 8. Acoustic pi -mode density kernels for mean (upper panel) and difference (lower panel) travel-time measurements between 
a point-pair separated by a distance of 25 Mm. The boundary of the spot, marked by the solid black line, is much smaller than 
the horizontal wavelength. The horizontal dot-dash line denotes the height at which observations are made in the quiet Sun 
and the symbols mark the measurement points. The pi-mode is seen to be much less affected by the spot than the / mode, 
as seen in the near symmetry of the kernels. The magnetic spot is indeed relatively weak and consequently, pi, whose modal 
energy is focused in the sub-surface layers, is much less sensitive to the field. 



